Fast Monte Carlo simulations and singularities in the probability distributions of 

non-equilibrium systems 
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A numerical technique is introduced that reduces exponentially the time required for Monte Carlo 
simulations of non-equilibrium systems. Results for the quasi-stationary probability distribution in 
two model systems are compared with the asymptotically exact theory in the limit of extremely small 
noise intensity. Singularities of the non-equilibrium distributions are revealed by the simulations. 

PACS numbers: 02.50.Ng, 02.70.Tt, 05.40.-a 



The understanding of fluctuations in systems away 
from thermal equilibrium is a problem of long standing in 
statistical physics |jy. Well known examples of physical 
phenomena in which non-equilibrium fluctuations play a 
particularly important role include e.g. switching of po- 
larization in lasers j^l], switching between different con- 
figurations in proteins the transition to instability in 
Josephson junctions P|, and chemical reactions 

In non-equilibrium systems, where symmetries of de- 
tailed balance are broken, no general methods exist for 
the calculation of even basic quantities like the prob- 
ability distribution. This is a case where numerical 
and asymptotic theoretical methods for investigating the 
probability distribution are of particular importance. 

Theoretical approaches, such as WKB-like or path- 
integral methods, are available in the limit of small noise 
intensity, ^ In particular the theory sug- 

gests that a solution to the problem of non-equilibrium 
fluctuations requires an understanding of the dynamics 
of deviations from the steady state I'l and an analysis 
of singularities in the non-equilibrium potential [Til Il2j . 
Some ideas about how to extend the existing (I? ^ 
limit) theory for still small but finite noise intensities 
have recently been suggested 0, 0, |23| ■ 

The main numerical technique used to verify theoret- 
ical predictions, and to analyse the behavior of the dy- 
namical system under study, is Monte Carlo simulation. 
The theory gives an asymptotically exact solution in the 
_D ^ limit. In contrast, D in the numerical simula- 
tions is necessarily finite. Typically, the time required for 
Monte Carlo simulations grows exponentially as D ^ 0. 
This meant that theoretical predictions of interesting sin- 
gular structures, and of the non-equilibrium probability 
distribution Il3l|, for long remained untested either ex- 
perimentally or by numerical simulation. Moreover there 
was no clear understanding of how the picture changes 
for small but still finite noise intensities. 

Approaches that have been tried to speed up the sim- 
ulations have focused mainly on finding optimal fluc- 
tuational paths and rates of transition between stable 
states of a system (e.g. efficient transition path sampling 
|l4j and dynamics importance sampling |15|. following 
the earlier suggestion of [3). In [13 the path sam- 



pling method was adapted for non-equilibrium systems. 
Based on the same idea, the umbrella sampling technique 
was suggested to estimate the probability of reaching any 
point in the phase space of an equilibrium system starting 
from a fixed initial state [14| |. A technique for improv- 
ing sampling in equilibrium systems by splitting up the 
probability packets was introduced in ;18]. So far, how- 
ever, no general algorithm has been suggested, able to 
give both the whole probability distribution and dynam- 
ical information like the optimal fluctuational paths for 
small noise intensities for non-equilibrium systems. 

In this Letter we introduce a numerical method that 
enables the time required for Monte Carlo simulations to 
be reduced by an exponentially large factor. It is applica- 
ble to generic two-dimensional non-equilibrium systems, 
does not require any a priori knowledge about the sys- 
tem apart from its dynamical equations of motion, and 
it allows the quasi-stationary probability distribution to 
be computed directly over the whole phase space. Using 
this method, we reveal for the first time singular behavior 
of the non-equilibrium distribution in numerical simula- 
tions, and we show that the results are in quantitative 
agreement with the asymptotic theory. 

The central idea is to perform the simulations in se- 
quential steps. We construct the quasi-stationary distri- 
bution, patching together intermediate results: we start 
from one of the steady states and gradually move away 
from it. Wc find that the time required for the simula- 
tions at each step is reduced by an exponentially large 
factor as compared to the standard technique: if the 
time necessary for a conventional Monte Carlo simula- 
tion technique is T, our modified method would require 



only time « NT exp 



-(AT-l). 



D , where TV is the num- 



ber of steps involved and is their separation in terms 
of the logarithm of the probability [2^ . 

We first explain the method on a very simple equi- 
librium stochastic system, and then we apply it to two 
much-studied non-equilibrium systems and compare the 
numerical results with theoretical predictions. To illus- 
trate the technique, we consider an overdamped Brown- 
ian particle moving in a bistable Duffing potential U{x) ~ 

-xy2 + xy4: 



x = -u'{x)+m, 



(1) 
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where ^(t) is zero- mean white Gaussian noise with inten- 
sity D and moments 

m) - 0, mm) = ^mt)- 

The probabihty distribution is completely defined by the 
potential U{x), and is of the Boltzmann form p(x) oc 
exp{—U{x)/D). As in the case of a non-equihbrium sys- 
tem (where the probability distribution is not defined by 
a potential) a standard Monte Carlo techni que can be 
used to deduce p{x). Numerical integration of the 
Langevin equation (^), assuming the system to be lo- 
cated initially at one of the potential minima Xm, gives 
the discrete probability distribution p{x), peaked at Xm- 
The potential can be deduced as $(x) oc —D\np{x). If 
the noise intensity is very small, the system fluctuates 
in a close vicinity of Xm and large deviations from it 
are extremely rare. Accordingly, the conventional Monte 
Carlo technique cannot be used to study the dynamics of 
optimal escape paths, or the properties of the probabil- 
ity distribution far from the potential minima: for small 
noise intensities the statistics required cannot in practice 
be collected within a realistic time. 

In order to overcome this problem, we start from the 
distribution already obtained We fix two proba- 

bility levels Pi and p/, lying well within the region where 
the numerical p is accurate, with pf < pi correspond- 
ing to two levels in the potential $i and and two 
coordinates Xi and x/, as shown in Fig. 1. We require 
the levels pi and pf to be fairly different, such that the 
corresponding Xi and Xf are sufficiently separated: the 
distance between them must exceed V Dh, where h is the 
integration time step used in the Monte Carlo simulation, 
and must also exceed the discretization step Ax in the 
discrete probability distribution. 

The next step of the simulation is started from the 
level $i (putting the system at its initial condi- 

tion) . If the system starts to evolve along a fluctuational 
trajectory (towards the boundary of attraction) we just 
follow its natural dynamics according to ^ and collect 
the statistics for building the probability distribution in 
the usual way. If the system starts with a relaxation tra- 
jectory (towards Xm), or when it crosses the boundary 
Xi due to relaxation some time later, we stop the simula- 
tion and reinject the system back to the initial state Xi. 
In this way we simulate the full dynamics of the system 
at higher levels of the potential ^{x) > ^i (in the re- 
gion of coordinate space x > Xi for this particular case) . 
Thus, in the subsequent simulation step we follow only 
those fluctuations that have already attained a certain 
level in the potential $i, without waiting for this expo- 
nentially slow event to happen. In this way, a new piece 
of the probability distribution is built with a time saving 
~ exp^i/D compared to a simulation starting from the 
potential minimum Xm- The computed new piece of the 
potential $2(2;) is shown as the upper curve in Fig. 1. 

The two pieces of the inferred potential (the original 
$1(2;) and the new $2(2;)) are then merged at Xf by a 
simple shift. Continuing this procedure, the probability 




FIG. 1; The first ($i(a;), lower curve) and second ("I>2(2;), 
upper curve) pieces of the inferred potential ^{x) for the sys- 
tem (0 with D — 0.005. The discontinuity in the gradient 
of 'I>2(a;) near Xi is an artefact due to a boundary effect in 
the calculation of the discrete probability distribution. To 
avoid this problem $i(x) and $2(3;) are merged at the point 
Xf and the initial part of $2(2:) is discarded. We normalize 
$1(2;) choosing $i(a;m) ~ 0, and each successive piece of $(a;) 
is normalized in order to match with the previous one at the 
point where they join. Inset: the inferred potential $(2:) for 
the system (0 with D = 0.005. The new technique (circles) is 
compared with standard Monte Carlo simulations (bold line) 
and with the Dufhng potential U (x) (thin line) . 

distribution and the corresponding potential can be built, 
step by step, for the whole region of interest. The inset 
in Fig. 1 shows the resultant potential, built from 13 
such pieces between the minimum at x^ = — 1 and the 
maximum at a: = 0. It coincides closely with the Dufhng 
potential U{x) itself. The potential ^{x) is thus inferred 
within a region of coordinate space that is inaccessible 
in a conventional simulation (shown as bold curve for 
comparison). We stress that no a priori knowledge of 
the dynamics has been used in the simulations, and that 
the method is robust to choice of parameters. 

In the case of a two dimensional system, the proce- 
dure remains essentially the same. The main difference 
is that, instead of identifying two points Xi and x/, we 
need to identify two closed lines of constant probability. 
One line is a boundary line for starting simulations from, 
and the other is a reference line for matching together 
different pieces of the probability distribution (see Fig. 2 
for clarification) ^25] . The crucial point of our technique 
is that, in starting the simulations from the boundary 
line, we must not perturb the natural dynamics of the 
system. This implies that we should consider the rein- 
jection location probability (RLP) along the boundary 
line corresponding to pi. Starting from the second step 
of the simulations, the system should be reinjected back 
according to the RLP after it relaxes across the bound- 
ary. We emphasize that the RLP is not the same as 
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the probability distribution p(x), which is constant on 
the boundary Hne. The RLP is an additional important 
measure which describes local discrete dynamics in the 
neighborhood of the boundary line. It is a distribution 
along the boundary of how often the system crosses it. 

In an equilibrium system, detailed balance provides a 
symmetry that can be used to reinject the system back 
at the boundary level, without any need to compute the 
RLP. For non-equilibrium systems, however, this proce- 
dure is inapplicable. The RLP should be considered sep- 
arately (and calculated explicitly) for the particular sys- 
tem being investigated. It can be obtained from an anal- 
ysis of the finite difference equation corresponding to the 
model. In the limit of small integration time step the 
probability to cross the boundary is proportional to the 
diffusion-related term in the finite difference equation. 
Then the RLP is simply proportional to the projection 
of the vector orthogonal to the boundary onto the coor- 
dinate affected by the noise ^. It can also be computed 
numerically. 

For non-equilibrium systems, the limit of small noise 
intensity is of particular importance. A sufficiently small 
D gives rise to the possibility of revealing the non- 
equilibrium potential 

<i>(x) = lim -£)lnp(x), 

directly through a numerical experiment. Observations 
of the predicted singular shape of Inp(x), and of its de- 
pendence on Z?, are thus of considerable interest. 








FIG. 2: The whole inferred $(a;,t) for the system Q for 
A = 0.1, n = 1, _D = 0.005. Two lines are the lines of 
constant probability found after the first step of simulations. 
The corresponding levels of probability were chosen as — 
3D and $/ = 5D. 

As a first example of a nonequilibrium system, consider 
the periodically-driven overdamped Duffing oscillator 

x = -U'{x)+Acosnt + ^{t). (2) 




-0.8 -0.6 -0.4 -0.2 jf. 



FIG. 3: A time section of the inferred ^{x,t = 4.1) for the 
system ^ with A — 0.1, i} — 1, and different noise intensities: 
D = 0.005 (diamonds); D = 0.01 (circles); and D = 0.02 
(triangles) . The theoretical predictions are shown by full lines 
for finite noise intensities, and by dashed line for D = 0. 
Inset; oscillations of ^{x,t) at the boundary of attraction for 
different noise intensities. 

We infer <I>(x,t) as —Dliip{x,t). This quantity corre- 
sponds to the theoretical "global minimum of the mod- 
ified action" in the Hamiltonian theory of large fluctua- 
tions |23| and, in the limit 13 ^ 0, it becomes the non- 
equilibrium potential. 

The complete $(a;, t), constructed from 12 such pieces, 
is shown in Fig. 2 and a time section of ^{x, t) calculated 
for different noise intensities together with the results 
of theoretical calculations (Hamiltonian theory including 
the prefactor) ^2d\ is shown in Fig. 3. The RLP in the 
simulations can be taken as constant if a small enough 
integration time step is used in the scheme. A small 
difference between the theory and the simulations results 
appears for larger noise intensities then the asymptotic 
theory starts to break down. 

As a second, more complicated, nonequilibrium exam- 
ple, consider the inverted Van-der-Pol oscillator 

X + 2ri{l - x'^)x + oj^x ^ ^{t) (3) 

Here, in order to be able to merge more easily the dif- 
ferent pieces of $(a;,y), we apply a coordinate transfor- 
mation from X and y ^ x to amplitude A and phase 
(j) {x = Acos{<j)), y — — Awo sin((/))). The probabihty 
p{x, y) can be then analyzed in the {A, (p) coordinate 
space. This makes the problem very similar to the pe- 
riodically driven Duffing oscillator: the only difference 
is the RLP which, in the case of the Van der Pol oscilla- 
tor, turns out to be strongly modulated. It is essential for 
this modulation to be taken into account when reinjecting 
the system back to the boundary of constant probability. 
Two sections of <^{x,y), obtained from the simulations 
for different parameters 77, are compared with the theory 
in Fig. 4. Again, the agreement between numerics and 
theory is excellent. 
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FIG. 4: A section {x = y) of the inferred '1>(^) for the system 
^ with ijjo = 1, noise intensity D — 0.01 and = 0.25 
(circles); and rj — 0.5 (diamonds). Theoretical predictions are 
shown in each case for D = (dashed curves) and D — 0.01 
(full curves). 

The non-equilibrium systems considered in this Letter 
share the same structure of singularities. Using the fast 
Monte Carlo simulations we reveal plateaus, the essen- 
tially flat regions in the probability distribution, which 
can be observed close to boundaries of attraction. They 
result from a purely dynamical effect that is not associ- 
ated with the flatness of any potential. We have shown 



that its origin is related to switching between different 
types of optimal fluctuational path, and it is a gen- 
eral feature of non-equilibrium systems with metastable 
states 20, The switching lines T^l are revealed as 
lines along which the "global minimum of the modified 
action" <i>(x) exhibits sharp bends - corresponding to 
the predicted line at which the non-equilibrium potential 
is non-differentiable. In the boundary region we found 
the oscillations of the probability distribution and their 
dependence on noise intensity (see the inset in Fig. 3) 
discussed in the recent publications 0, 0, l^^- Using 
the simulations we demonstrated noise induced shift of 
the singularities and the optimal escape path, which has 
stimulated a new step in the development of the theory 

We emphasize that the singularities can be confidently 
observed only in the limit of extremely small noise inten- 
sity, and therefore that the use of our new technique is 
crucial in that it reduces by an exponentially large factor 
the time required for Monte Carlo simulations. In ad- 
dition to being fast, it preserves dynamical information, 
can be modified to analyse optimal fluctuational paths, 
is applicable to the energy-optimal control problem , 
and can be further extended to encompass higher dimen- 
sional systems and maps. 
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